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Abstract 

Careful studies of liquid water in its supercooled region have led to many insights into the struc- 
ture and behavior of water in all temperature regions. While bulk water freezes at its homogeneous 
nucleation temperature of approximately 235 K, for protein hydration water, the binding of water 
molecules to the protein avoids crystallization. Here we study the slow dynamics of water at low 
temperatures by dielectric relaxation experiments on lysozyme hydration water over an extremely 
broad range of frequencies and temperatures. We probe the temperature dependence of the re- 
laxation time due to water protons in order to gain insight into the structure and connectivity of 
the hydrogen bond network of water molecules adsorbed on the protein surface. We observe two 
dynamic crossovers: one at about 252 K, assigned to a change of the diffusive regime of water 
protons along the hydrogen bond network; the other, at about 181 K assigned to a structural 
re-arrangement of the water network. We support this interpretation by comparing the experi- 
mental results with Monte Carlo simulations and mean-field calculations on a cell model of water. 
The model allows an interpretation in terms of thermodynamic features of water, predicting the 
presence of two specific heat maxima at ambient pressure. The first is due to fluctuations in the 
hydrogen bond formation, and the second, at lower temperature, due to cooperative reordering of 
the hydrogen bond network. 

PACS numbers: 
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Among liquids, water is undoubtedly the most studied, as it is the medium of life and 
chemistry on our planet, and its puzzling properties are a challenge for science. Many 
aspects of the dynamical behavior of water, both as a bulk solvent or adsorbed on the 
surface of biological macromolecules, can be tackled by dielectric relaxation spectroscopy. 
This powerful technique measures the interaction of a sample with a time-dependent small- 
amplitude electric field [l| and is one of the few techniques that can span a very large 
frequency range, allowing the observation of dynamical processes over many time scales. 
The field-induced polarization can be investigated either in the time domain or as a function 
of field frequency, u. In the latter case the response function includes all the contributions 
to polarization that depend on frequency, and thus reflects the dynamics of dipole and the 
displacement of charges in the sample. Recent work has focused on dielectric relaxation 
experiments on hydrated proteins, aimed at the study of the dynamics of water molecules 
together with the protein amino acids . At variance with these experiments, we will 
report on the dielectric relaxation of water protons. In particular, we discuss the temperature 
dependence of proton dynamics along chains of hydrogen bonded water molecules at the 
interface with a globular protein. The rationale behind this approach is that measurements 
of the proton mobility are closely linked to the dynamics and connectivity of the water 
molecules in the protein hydration shell {4]. The dielectric relaxation of water protons is a 
sensible probe for the structure of the underlying hydrogen bond network, and therefore of 
particular interest in the present context. Moreover, protein tumbling can be neglected if 
sample hydration level is kept below 0.4 g H 2 0/g dry protein, and results can be interpreted 
in the framework of well established theoretical models 

A dielectric spectrum of the hydrated lysozyme sample at 215 K is shown in Fig. [IJ 
The high frequency dispersion can be compared to that observed in a recent dielectric 
study of the same hydrated protein 2I, and we find an identical temperature dependence 
(data not shown). The assignment of this relaxation process is still controversial {2I, 8|. In 
addition, results obtained with other experimental techniques applied to hydrated proteins 
have raised an intense debate on the nature and temperature dependence of this dynamical 
process (9), lo| . 

Here we focus on a secondary relaxation process, previously tentatively ascribed to motion 
of protein side chains or to a cross-term due to the relaxation of water and protein dipoles . 
We have shown [6j that this secondary relaxation can be resolved into two contributions 
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(Fig. [I]): one effectively due to the relaxation of protein side chains, and the other to the 
motion of water protons alongthe hydrogen bond network of water molecules in the first 



hydration shell of the protein 



111- 



Figure [2(a) shows the temperature dependence of the water proton relaxation time r for 
0.3 g H 2 0/g dry protein over the temperature range 170 K to 300 K. This relaxation time 
has a very strong hydration dependence 12] , and in lysozyme samples with hydration lower 
than 0.28 g H^O/g dry protein it can be measured only above 230 K. 

At high T, t shows a non-Arrhenius temperature dependence, well described by the 
Vogel-Fulcher-Tamman (VFT) equation r(T) = r exp [Bt/R(T — To)], where r , Bt and 
T are fitting parameters and R is the gas constant. At about 252 K, a crossover is observed 
(Fig. [2(a)), from a non-Arrhenius behavior to a second non-Arrhenius behavior. 

It should be noted that at about the same temperature, namely 254 K, a dynamical 
transition was found for the same hydrated protein, and described in terms of random walk 
of protons along the water network changing from a sub-diffusive regime at low temperature 
to a diffusive regime above T = 254 K [7|. 

At about 181 K we observe a second crossover (Fig. [2(a)), as the temperature dependence 
of the water proton relaxation time changes from a VFT above 181 K, to Arrhenius, r(T) = 
To exp (A/RT) below 181 K, where To is a characteristic relaxation time and A is a constant 
activation energy. This might be due to a slowing down of the proton mobility with a less 
marked temperature dependence, possibly due to a structural re-arrangement of the water 
network. 

Neutron scattering experiments on the same hydrated protein revealed two dynamical 
transitions, one at 220 K and the other at 150 K 13J. These transitions have been attributed 



respectively to the rotational motion of interfacial water, and to proton dynamics on a local 
(few A) scale. In particular, the low temperature transition (at 150 K) was claimed to be 
due to a sudden increase of the configurational entropy of the system, linked to a significant 
excursion of the hydrogen bond length [lsj]. Deep inelastic neutron scattering experiments 
on the same hydrated protein question this interpretation, as they do no reveal any change of 
the hydrogen bond length and of the proton potential up to 290 K ll|, [l4j. 2 H-NMR studies 



on hydrated proteins at a comparable hydration level as that investigated in the present 
report have nevertheless indicated two dynamical transitions at about 250 K and at 180 K 
due to protons in the protein hydration shell [h]]. While the low temperature transition 
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has been interpreted in terms of local motion of protons in a rather rigid hydrogen bond 
network, consistent with the present work, the high temperature transition is not discussed 
there in the same framework |l0|. Our goal here is to give a consistent thermodynamic 
interpretation of both crossovers. 

To gain insight into the behavior of the relaxation dynamics we perform Monte Carlo 
(MC) simulations and mean-field calculations of a cell model for water. Since the protein 
powder is, on average, hydrated only by a single layer of molecules, and lysozyme is a 
compact globular protein which does not undergo any configurational transformation in our 
experimental conditions, we approximate that the only effect of the water-protein interaction 
is to reduce the hydration water to a d = 2-dimensional system with the water coordination 

n 

number fixed to four [15[. As a consequence each water molecule can form at most four 
hydrogen bonds. Since we focus on the hydrogen bond dynamics of water molecules, for 
our scope is not relevant if a hydrogen bond is formed with the protein or between water 
molecules, but it is relevant that water molecules can form a network of hydrogen bonds. 

We consider N = 10 4 molecules in a d = 2-dimensional volume V, with periodic boundary 
conditions, divided into N equivalent cells each with one molecule % G [1, N] and with volume 
V/N larger than a hard-core volume yn due to van der Waals repulsion. The interaction 



Hamiltonian in the liquid phase is 16|, [l7| 



= Uo (r) - JJ2 5 ° m - 3 ° E 5 — • W 

(ij) (k,l)i 

The first term Uo(r) denotes the sum of all the isotropic interactions (such as the van der 
Waals interaction) between molecules at distance r = (V/N) 1 / d , represented by a Lennard- 
Jones potential with attractive energy e plus a hard-core at distance r$ = (vq) 1 ^. Compar- 
ison of our simulations results with the experiments presented here allow us to set e = 5.8 
kJ/mol, consistent with the value 5.5 kJ/mol of the estimate of the van der Waals attraction 



based on isoelectronic molecules at optimal separation 18), and r = 2.9 A consistent with 



3, 
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the expected value of the van der Waals radius 

The second term (with J = 2.9 kJ/mol, 5 a ^ = 1 if a = b and 5 a ^ = otherwise, and 
(i, j) denoting that % and j are nearest-neighbors) accounts for the directional contribution 
to the hydrogen bond energy, where cr^- = 1, . . . , q is a (Potts) variable representing the 
orientational state of the hydrogen (or the lone e~) of molecule % facing the lone e~ (or the 
hydrogen, respectively) of the molecules j. For the sake of simplicity we do not distinguish 



between hydrogen and lone e~, associating to each molecule four equivalent bond indices 
(Jij. Assuming that a hydrogen bond is formed when the angular deviation with respect to 
the linear bond of the OOH angle is less than 30°, we get q = 180°/30° = 6. The third 
term (with J a = 0.29 kJ/mol, and (k,l)i indicating each of the six different pairs of the 
four bond indices of molecule i) represents an interaction accounting for the hydrogen bond 
cooperativity and giving rise to the 0-0-0 correlation [21], locally driving the molecules 
toward an ordered (tetrahedral in the bulk) configuration with lower energy. By defining 
the energy per hydrogen bond (between cr^ and <jji) as the sum of the interactions in which 
two bonded molecules {i and j) are participating, we get Ebb = e + J + mJ a /2, where 
m = 0, . . . , 6 is the number of cooperative interactions in which that bond variables (cr^- and 
<jji) are participating. With our choice of parameters the values of ranges between 8.70 
and 9.6 kJ/mol depending on m. However, a definition of .Ehb based on a cluster of 5 or 8 
bonded molecules in 3D increases this range up to 17 or 18 kJ/mol, respectively. Therefore, 



E'hb depends on the environment (the value of m and the number of molecules 



Donded in 
and has 



a cluster), as observed in computer simulation of the crystalline phases of ice [22 
values within the range 6.3 H - 23.3 kJ/mol [2J], proposed on the base of experiments. 

The total volume of the system V = Vmc + -Whb^hr can change to accommodate the 
volume expansion v hb due to hydrogen bond formation [25| , where Vmc ^ is a dynamical 
variable allowed to fluctuate in the simulations, and A^b = J2<i j> ^aava * s ^ ne total number 
of hydrogen bonds. We adopt Whb = 0.5wq, correspondin g to a maximum hydrogen bond 
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distance of about 3.3 A, consistent with other authors 

We perform MC simulations at constant N, T and pressure P (see Methods). At low T 
and high P we find a liquid-liquid first order phase transition, with negative slope in the 



P-T plane, ending in a critical point at P c = 0.18 ± 0.01 GPa and T c = 171 ± 1 K [17|, |27j. 

We compute the autocorrelation function Cjv/(i) (see Methods) of the quantity that gives 
the orientational order of the four bond indices of molecule i. In Fig. E|b) we show the 
calculation of the MC relaxation time, defined as the time at which Cjvr(t) = 1/e, as a 
function of 1/T for pressure P = 0.1 MPa ~ 1 atm. We find two crossovers: a non- 
Arrhenius to non-Arrhenius crossover at T w 252 K and a non-Arrhenius to Arrhenius 
crossover at T « 181 K. 

From general considerations we know that the structure of a liquid must influence its 
dynamics. The Adam-Gibbs theory 28J] relates a maximum in isobaric specific heat Cp as 
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a function of T with a crossover in relaxation times. This can be intuitively understood if 
we consider that Cp = T(dS/dT)p corresponds to the variation of entropy S at constant 
P, (dS/dT)p. When Cp is large, S has a large decrease with decreasing T, i.e. the amount 
of order increases and the liquid becomes more structured. Since the dependence of the 
relaxation time r on T is stronger when the energy change needed to modify the structure 
is larger, and the change of energy is proportional to Cp, a maximum variation in the 
structure (Cp 1 *™ ) implies a change (crossover) in the way r depends on T. Recent computer 



simulations have exploited this relation 23. l30l|. 



We calculate the isobaric specific heat Cp = (dH/dT) p , where H = (E) + P(V) is the 
enthalpy, and (•) denotes the thermodynamic average. Figure [3(a) shows Cp as a function 
of T for P = 0.1 MPa. Two maxima in Cp are clearly visible, at variance with previous 
results 31|. As the liquid-liquid critical point is approached in P, the narrow maximum 
at lower T diverges, while the relatively broader maximum does not, and the two maxima 
merge. 

To understand the origin of the two Cp maxima, we write the enthalpy as the sum of two 
terms 

H = H HB + H Coop (2) 

where H HB = (-JN HB + P(V MC + N K bVhb)) and F Coop = H - H RB . Hence, we consider 
Cp = Cp B + Cp° op where we define the component C™ = [dH RB /dT) and the coop- 
erative component Cp° op = (dH Coop /dT) p [Fig. [3(a)]. Cp B is responsible for the broad 
maximum at higher T. Cp B captures the enthalpy fluctuations due to the hydrogen bond 
formation given by the terms proportional to the hydrogen bond number A^hb- To show this 
relation, we calculate the locus of maximum fluctuation of (Nub), related to the maximum 
of |d(A^HB)/dT|p [Fig. [2(b)], and find that the temperatures of these maxima correlate very 
well with the locus of maxima of Cp B . 

We find in Fig. [3j(a) that the maximum of Cp at lower T is given by the maximum of 
C p oop . To confirm that Cp° op corresponds to the enthalpy fluctuations due to the cooperative 
term in Eq. [TJ proportional to J a , we calculate |d(iV Co0 p)/dT|p, where (N Coop ) is the average 
number of molecules with complete local order of their bond indices. (In the range of T of 
interest here the contribution to H of the Uo term is negligible.) We find that the locus of 
maxima of |diVcoop/dT|p [Fig. [31(b)] overlaps with the locus of maxima 

ofC ,Coo P The 

same 

qualitative behavior is predicted on the base of mean-field (MF) calculations for the model 



[Fig. He)]. 

The Arrhenius behavior of r at very low T is, therefore, the consequence of the saturation 
of the hydrogen bond network (|d(iVHB)/dT|p ~ 0). In both experiments and simulations, 
the constant Arrhenius activation energy A for structural rearrangements (Fig. [2]) is consis- 
tent with the average necessary to break a hydrogen bond at low T. In the high T limit 
instead, the non-Arrhenius behavior of r is a consequence of the dilution and disorder of 
the hydrogen bond network and the activation energy increases for decreasing T as does the 
value .Ehb averaged over all the molecules. At T just below the maximum of |d(A r H B)/dT| P , 
the rate of change of (-Ehb) decreases as does |d(iVHB)/dT|p, but at lower T it increases again 
as a consequence of the increase of |d(A^ Co op)/dT|p, giving rise to another non-Arrhenius 
behavior of r at intermediate T, up to reaching the maximal fluctuation, below which the 
hydrogen bond network saturates and the Arrhenius behavior takes place. 

The qualitative agreement with mean field results (where the dimensionality d and the 
nearest-neighbor distance play no role and where the microscopic changes of the water 
structure induced by the protein matrix are reflected only in the maximum number of 
nearest-neighbors, here fixed to four) shows that we can expect the same behavior also for 
the model in d = 3. Therefore, the model predicts that what is observed in the experiments 
could be a property of water (including bulk water) and that the presence of the protein 
matrix is relevant because it induces a modification in the structure of water inhibiting 
crystallization, and makes low T observation possible. 

Our results indicate that the thermodynamic features of water induce two dynamical 
crossovers. The fluctuations of the hydrogen bond formation give rise to a broad maximum 
in Cp, which corresponds to the non-Arrhenius to non-Arrhenius crossover, at higher T. 
The cooperative ordering (restructuring and saturation) of the hydrogen bond network, 
corresponding to the sharp maximum of Cp at lower T, gives rise to the non-Arrhenius 
to Arrhenius crossover. Our preliminary results show that the effect of the hydration level 
and the pressure is to increase the cooperative behavior of the hydrogen bonds, reducing 
the range of T where the intermediate non-Arrhenius behavior is observed and reducing the 
separation in T between the two Cp maxima. At high hydration level in the experiments, 
or at high P in the model, the intermediate dynamic regime disappears and the two Cp 
maxima merge. In the model the two maxima of Cp approach each other by increasing P 
below the liquid-liquid critical P, while progressively separate at P higher then the liquid- 
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liquid critical point. 



METHODS 

-Experiments: Measurements were performed using an Alpha Analyzer apparatus (Novo- 
control) in the temperature range from 170 K to 300 K at constant temperature steps and 
in the frequency range from 10~ 2 Hz to 10 7 Hz. Powder samples of crystallized and highly 
purified lysozyme powder (from chicken egg white) were purchased from Sigma-Aldrich. 
Sample pH and hydration h were adjusted to a final value of 7 and 0.3 g H 2 0/g dry pro- 
tein, respectively, as previously described [6]. Data analysis, described elsewhere p, was 
performed fitting the measured complex spectra, taking into account sample permittivity 
(expressed as a combination of Havriliak-Negami functions and a conductivity term) and 
electrode polarization effects 

M 



7=1 I s - 



where cq is the conductivity, is the high frequency limit of the permittivity, M is the 
number of relaxation processes, Aej and Tj are the dielectric strength and the relaxation 
time for the jth contribution, respectively, and A and d characterize the polarizability term. 
-Simulations: The MC dynamics consists in updating the variables cr^ by means of 



the Wolff algorithm 32J and updating the volume V in accordance with the acceptance 
probability min (l,exp [-0 (AE + PAV - NRTln(V f /Vi)))), where (3 = (RT)~\ VJ and V f 
are respectively the initial and final values of the volume, AV = Vf — VJ, AE = Ef — Ei 
and Ei and Ef are respectively the initial and final values of the energy calculated as the 



the right hand side of Eq. [T] [27]. 



The autocorrelation function for the quantity Mj = | ^ . cr^ , that gives the orientational 
order of the four bond indices of molecule i, is defined as 

r m - 1 ^ (M t (t + t)M t (t ))-(M t ) 2 

CM{t) = N^ ( Ml (t y) - (Mi) 2 • (4) 

To rescale the MC results on the experimental results we apply a linear rescaling on T and 
P: for a temperature expressed in MC internal units T[e/R], the corresponding temperature 
in K is T[K] = (134.5 + 700 T[e/R]) K, corresponding to e = 5.8 kJ/mol, and for the 
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pressure P[GPa] = (7.56 x 1(T 4 + 0.257P[e/v }) GPa, corresponding to r = 2.9 A. For the 
correlation time r we apply a logarithmic rescaling m(r[s]) = —31.3 + 1.741nr[MC Steps], 
that is consistent with the observation that the conversion from MC Steps to real time is 
T-dependent, with a MC Step corresponding to a larger time at lower T. 

-Mean-field (MF) calculations: Each cell is assumed to have constant volume t>o, an d 
is described by occupancy variable = 0, 1. The system is assumed homogenous, with 
average occupancy n = Yl,i n i/N- With a cluster expansion we approximate the probability, 
p a (h a , T, P), of forming a hydrogen bond between neighboring molecules, expressed in terms 
of P, T and the hydrogen bond density field h a . MF expressions for the molecular entropy so 
of the variables n^, and s a of the variables <Jij, are derived. In equilibrium at constant T, P, 
and chemical potential fi, the system will tend to minimize the Gibbs free energy. Hence the 
molecular Gibbs free energy, g(h a ,T,P) = —2en—2(J~PvHB)'npa—6JaPa+P'Vo—T(s + s a .), 
is minimized with respect to n and h a at each (T, P), yielding equilibrium values n eq (T, P) 
and h eq (T,P). Substitution of n eq and h eq into the MF expression for the volume V = 
Nv + 2Nn 2 p a gives the full equation of state of the system. 

In MF we rescale T with the experimental values as T[K] = (102.5 + 850T[€mf/-R]) K, 
corresponding to 6mf = 706.7 kJ/mol. 
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T=215 K 

side chain relaxation 

proton relaxation 
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FIG. 1: Dielectric spectrum of hydrated lysozyme sample at 215 K. The sample response function 
is expressed in terms of the complex permittivity e^(a») = e'(uj) — ie"(uj). This figure shows the real 
(e', • left axis) and imaginary (e", A right axis) components of e* m as a function of angular frequency 
lo. The solid lines trough the symbols are the result of the fitting procedure in the complex plane 
O]. The three peaks related to the main relaxation, to the water proton relaxation, and to protein 
side chains relaxation, contributing to the measured imaginary component, are shown. From this 
figure the contribution due to sample conductivity and the term due to electrode polarization (see 
Eq. [3]) are omitted for clarity. These terms contribute to the low frequency tail of the measured 
complex permittivity. 
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FIG. 2: Two crossovers in the water proton relaxation time r: a non-Arrhenius to non-Arrhenius 
crossover at T ~ 252 K and a non-Arrhenius to Arrhenius crossover at T ~ 181 K. (a) Experimental 
temperature dependence of r (•), as a function of inverse temperature (1/T) at ambient pressure. 
Solid line is a fit to a VFT equation r = To exp[Bx/R(T — Tq)], with fitting parameters To = 
7.8 x 1(T 12 s, B T = 9.4 kJ/mol, T = 180 K, and describes t in the high T region. Dotted line 
indicates the VFT temperature dependence in the intermediate T region, with fitting parameters 
To = 6.5 x 1CP 8 s, Bt = 6.2 kJ/mol, To = 140 K. Dashed line is a fit with Arrhenius equation 
t = Toex.p[A/RT] with fitting parameters to = 1.1 x 10~ 7 s, A = 25.2 kJ/mol, and describes the 
T dependence of r in the low T region. Notice that the two limiting behaviors, namely that at 
high T and that at low T, intersects at about 232 K. (b) Relaxation time (•) from MC calculations 
as a function of 1/T, for P = 0.1 MPa. The solid line is a fit to a VFT, with fitting parameters 
t = 1.61 x 10~ 8 s, B T = 5.2 kJ/mol, T = 181.2 K; the dotted line is also a fit to a VFT with 
parameters To = 7.5 x 10 -10 s, Bt = 15.9 kJ/mol, To = 95.2 K. The fitting parameters for the 
Arrhenius law (dashed line) at the lowest T are To = 3.3 x 1CT 4 , A = 13.7 kJ/mol. In both panels 
the Arrhenius activation energy is consistent with the energy to break one hydrogen bond in the 
experiments or the model, respectively. 
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FIG. 3: The isobaric specific heat Cp from MC simulations for P = 0.1 MPa. (a) Cp (•) can 
be decomposed into the component Cp° op (o) due to the hydrogen bond cooperativity, and the 
component Cp B (□) due to fluctuations of hydrogen bond formation, (b) The variation with 
temperature of the average number of cooperative bonds (iVcoop) and of the average number of 
hydrogen bonds (-/Vhb) for P = 0.1 MPa. (|d(A?c oop )/dT|)p and (|d(iV HB }/dT|)p show maxima 
where Cp has maxima, (c) Comparison of Cp from MF calculations for the case without cooperative 
interaction (J CT = 0, labeled as Cp B ) and for the case with cooperative interaction (J a /e = 0.05, 
labeled as c£ oop + Cp 16 ). The low-T maximum is present only in the second case, hence is due to 
the cooperativity. Calculations are at P = 0. 
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